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METHOD OF USING IMATRIX RANK REDUCTION TO REMOVE RANDOM NOISE FROM 

SEISMIC DATA 

FIELD OF THE INVENTION 

The present invention relates generally to processing seismic data and particularly to 
redudng noise in seismic data using a variety of 3D eigen filtering techniques based on matrix 
rank reducKon in the frequency domain. 

BACKGROUND OF THE INVENTION 

Stismic data can be used to interpret or to infer sub-surface geology, making it useful for 
the location, kientification and exptoitation of petroleum and minerals. However since seismic 
traces are often contaminated by random noise, seismic data sets typically undergo a series of 
conventbnaf statistical processes (known as "seismic processing") before the data can be so 
used. It is advantageous to remove noise at an earty stage in processing, since this improves 
the ability to perform all subsequent processing woric. Conventional seismic processing 
disadvantageousty degrades or distorts signal prior to removing sufficient noise for some 
purposes. Conventional seismic processing methods are based on functions and their 
transforms that it is often useful to think of as occupying two domains. These domains have 
been referred to as the function and transform domains, but more commonly they are known as 
the time and frequency domains. Operattons performed in one domain have corresponding 
operattons in tiie other. For example, the convolutk)n operation in the firm domain becomes a 
multiplication operation in tifie frequency domain and the reverse is also true, permitting users to 
move easily between domains so that operations can be performed where they are easiest. A 
seismic audk) signal being sampled at 8 Hz means that at each successive eightii of a second a 
measurement of the intensity of the signal is taken. The Fourier tran^orw decomposes or 
separates a wavefomn or function into sinusoids of different frequency, which sinusoids sum to 
the original waveform. It kientifies or distinguishes tite different frequency sinusokis and their 
respective anrv)litudes. The Discrete Fourier Transform (OFT) is useful because a digital 
computer works only with discrete data, numerical computation of the Fourier transform of f(t) 
requires discrete sample values of f(^. which are called 4. In addition, a computer can compute 
tiie transform F(s) only at discrete values of s. in other words it can only provide discrete 
samples of the transform. It is well understood tiiat a Discrete Fourier Transform (*'DPT) is a 
Fourier transform calculated for a wavelet over a finite interval so that values are given only for 
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the fundamental frequency (redpnxal of the inten^ai) and its harmonics. And a Irace** is a time 
domain record of seismic data from a single eiectronrtagnelic channel. 

Within the seismic industry sections of data are comprised of several "grids- that form 
cubes of data when used to correlate Cartesian (x-y surfece) coordinates to either a time 
domain or frequency domain that covers a finite range of time or frequency respectively. In a 
Time Cube; a trace Is a temporal record of date from a single seismic tihamel the sensor for 
vMidt\ is located on a point of intersection defined according to a given scale of abscissas (i.e. 
along the X-axis) and ordinates (i.e. along the Y-axis) - each of which points of intersecnon are 
identified by a co-ordinate pair and all of which points of intersection form a planar rectangular 
"grW". The 3"* dimensiwi being time, a "Irace" forms as seismic energy (typically in the fomt of 
acoustic signal amplitude) sample recordings are captured from the sensor location identtfied by 
the sutsject co-ordinate pair, over a finite range of finne (e.g. 2 seconds) at defined (e.g. 10 
millisecond) intervals. The scope and depth of the data comprising a Time Cube is related to the 
number of elements (Cartesian co-ordinate pairs), the time, and the sampling rate through that 
time. For example, in each 25 x 25 **grid'* there will be 625 traces forming a time record of 
signals (comprising both genuine reflected energy and noise) recorded in relation to the subject 
section (i.e. 1 date element assodated with each co-ordinate pair « 625 date elements per 
sample period). With each trace having samples captured every 10 milliseconds over 2 
seconds, this Time Cube would contein 125,000 date elemente. 

Prior Art Figure 1 illustrates a typical gedogical exploration setup for acquiring seismic 
data. Positioned at or just below the earth^s surface 110, an energy source 120 (typically 
explosive high-energy or vibrational low energy source having an effective length of only 20 to 
40 ms) generates at least one sound wave or signal (the wavefront resulting from which is 
typically recorded for approx. 3 seconds) having sufficient energy to follow path 130 down into 
the earth a suiteble distence to reflect at the interface of any changes in geology (commonly 
known as ^'events" or "rsftectors") 140, the reflected energy traveling beck to the surtece via 
path 150 is stmulteneously recorded by receivers 160 (commonly geophones positioned as an 
anay for 3D. or a line for 2D exploration). In marine seismic the sound wave is generated just 
below the surtece of the water and the reflected energy detected by hydrophones. For each 
such sound wave generation, or "shoT, the reflected signal returning to the surface via path 150 
creates a Irace'. which is a single recording at a receiver 160. Traces are detected and 
recorded in the form of a time series of sample measurements of particle velodty (for tend date) 
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or pressure change (for marine data). Many shots are taken to generate a seismic data set, 
often resulting in hundreds of nf^illions of traoes that may be stadced of summed in a variety of 
ways. When high energy impulsive seismic sources are used, the level of the detected true 
earth response seismic signal is usually greater than the ambient noise. However, when low 
energy surface seismic sources are used, the ambient noise can be at a level greater than the 
true earth response seismic signal. For this reason, seismk>trace recordings are often made 
involving the repeated initiation of a tow energy surface seismic source at about the same 
origination point, theret>y producing a sequence of sefsmio-trace data based on seismic wave 
reflections and/or refractions that have traveled over about the same path and therefore have 
approximately the same travel times. The process of adding such seismio-traoe data together 
for improving the signal-to^oise ratio of the composite seismic-trace recording is Icnown as 
"Vertical compositing" or Vertical stacking." it should be distinguished from -horizontal stacking " 
a process applied to a sequence of seismic-trace data based on seismic wave neflections from 
approximately the same subsurface point (referred to as the "common-depth point;* or "CMP") 
but which has been generated and recorded at different surface locations. Horizontal stacking of 
CMP seismiC'traoe data requires that time corrections (called "nonnal moveout" or -NMO; 
conrections) be applied before the traces are summed together, since travel times from seismic 
source to detector are unequal for each trace in the sequence, (t can be assumed that the true 
eartti response seismic signal embedded in each trace is coherent and in phase (con'eiated) 
and that ttie noise is random and incoherent (uncorrelated) witt> zero mean value. In general, 
the (^jective of vertical stacking is to maximtze Qie signal-to-noise ratio of the resultant 
recording. Reflectors that are not Haf are said to "dip" or slope. 

The use of a low energy vibrator can be more economical than Vt\e use of dynamite. 
Furthermore, as compared to the use of a high^energy impuls'ive seismic source, such as 
dynamite, the frequency of the seismic waves generated by a vibrator can be selected by 
controlling the frequency of tire pilot signal to 9\b power source, such as a hydraulic motor, 
which drives the vibrator. More particulariy, the frequency of tfie pitot signal to the vibrator power 
source can be varied, that is, "swept," for obtaining seismic-trace data at different frequencies. A 
low energy seismic wave, such as generated by a vibrator, can be used effectively for seismic 
prospecting if the frequency of the vibrator "chirp** signal which generates the seismic wave is 
swept according to a known pilot signal and the detected seismic wave reflections and/or 
refractions are cross-correlated with the pitot signal in order to produce seismic-trace recordings 
similar to those which woukl have been produced with a high energy impulsive seismic source. 
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Typically, the pilot signal is a swept frequency sine wave that causes the vibrator power source 
to drive the \ftbrator for coupling a swept sine wave "chirp" signal into the earth. The swept 
frequency operation yields seismi'c^trace data that enables different earth responses to be 
analyzed, prodding a basis on which to define the stnjcture. such as the depth and thldcness, of 
the subsurface fbmnations. It is a problem that recorded seismic-trace data always includes 
some background (ambient) noise in addition to the detected seismic waves reflected and/or 
refracted from the subsurface fomiations (referred to as the "tme earth response"). Noise can 
be classified as "stationary" and "non-stationary*, both of which can be random. Stationary noise 
is random n(»se suc^ as atmospheric electromagnetic disturbances that are statistically time- 
invariant over the period of acquisition of seismic-trace data for producing a recording. Non- 
stationary nc»se is random and often occurs as bursts or spikes generally caused by wind» 
traffic, recorder electrical noise, et cetera, which are statistically time-variant over the period of 
acquisition of seismkHraoe data for producing a recording and exhibits relatively large 
excursions in amplitude. In connection with swept frequency operation of low energy \ribrator 
seismic prospecting, it is common practice to vertically stack, or sum. the seismio^traoe data 
from a series of initiations, that is. sequential swept frequency operations, to produce a 
composite seismic-trace recording for the purpose of improving the srgnal4o-noise ratio of the 
seismic-trace data. Unfortunately, the commonly used technique of vertically slacking trace data 
is Inadequate in ttie presence of non-stationary noise ttiat appears during such seismic 
prospecting. 

Seisn\k; data is acquired in two principal geometries: 2-D and 3-0. In 2-0 acquisition, 
shots and recovers are positioned atong a (not necessarily straight) surface line. In 3-D 
acquisition, shots and receivers are positioned over a 2-D surface area. Seismic data related to 
3D geologic volumes necessarily includes random noise that may be isc^ated from Vdb signal 
data to different degrees by different conventional techniques, including an eigenimage filtering 
technique that is 20 in nature and disadvantageously does not account for additional available 
information respecting the formation. 

For 2-D acquisition the main product of seismic processing is a 2-0 stacked "section" 
(illustrated in Prior Art Figure 2). one of the dimensions representing horizontal position atong 
acquisition line 210, and the otiier dimension representing time 220. For 3-0 acquisition, 
seismic processing resulting in a 3-0 stacked section (illustrated in Pnor Ait Figure 3), two of the 



SEP. 12.2003 1:1?PM 60WUMGS 1403 263 9193 " " N0.7531 --P. 



dimensions representing edges 310 and 320 of the acquisition surfaoe area, and the o^er 
dimension representing time 330. 



Known seismic processing steps include commonHfnidpoint (CMP) stacking, where 
traces are collected into groups having roughly the same midpointe between the locations of the 
shot by whi<^ they were gyrated and the receiver at whidi ttiey were detected. For each 
recorded time sample, the magnitudes or values of every trace in the group are summed 
together, produdng a single "stacked" trace for each group. Such stacking commonly reduces 
the amount of data that must be processed by a factor of 10 to 100. 

Geological interpretatton is easiest and most successftjl on seisnrtic sections having low 
levels of noise, and thus one of the objects of seismic processing is to remove as much noise as 
possible. Noise can be broadly categorised as random, coherent, or monochromatic. Random 
noise may be defined as noise that is uncorreiated between traces and spectrally broad band. 
Some of the causes of random noise are the effecte of wind and other disruptions on the 
seismk: receiver and cable, poor penetration of seismic energy through the earth (part'cuiariy in 
the near surface beneath the shot or receiver), and numerous natural and man-made seismic 
energy sources apart from the intended one. The most conmon stage to carry out random- 
noise removal is after CMP stacking. A number of mettKXis have been devek>ped to do this, 
including: f4< transform (March and Bailey, 1983), f-x predlctkm (Canales, 1984; Soubaras, 
1994), Karhunen-Loeve transform (Jones and Levy, 1987; Al-Yat^, 1991). etgenimage 
(Ulrych, Saochi, and Freire, 1999). spectral matnx filtering (Gounon, Marse, and Goncalves. 
1998), and Radon transfonn (Russefl, Hampson, Chun, 1990(a) and 1990(b)). 



The foregoing methods work on 2«0 data sets, but can be adapted fbr 3-D stacked 
secSons by slk^ng the data volume along one of the spatial dimensions, Ntering each of these 
slices separately as if it were a 2-D sectton, and then recomposlng ttie 3-D volume. This can 
0ien be repeated in the opposite spatial direction. Such methods are not optimum in that they 
fail to fully expkxt the large amount of data available within a short radius of each spafaal point of 
the 3-D volume. For this reason, 'true S-D*" methods have been devetoped that work in both 
spatial dimensions simultaneously, including: f-xy prediction (Chase, 1992; Soubaras. 2000) 
and f4d( transform (Peardon and Bacon, 1992). 
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f^ndom noise removal before CMP stacking is less common. There are. however, at 
ieaet two advantages to removing random noise as early in the processing stream as possible. 
First it improves the performance of subsequent processes, notably deconvolution, statics 
conrecUon, and velocity analysis. Second, it has the potential to be more effective since more 
data is available before stacking, providing better statistk:al redundancy. At the same time, extra 
data means that random noise removal before CMP stacking requires more computation. Ndse 
removal before statics or deconvolution faces the problem of "surfece-consistent effe(^s^ 
meaning effects that are constant within each shot and receiver, but that may change radically 
even between adjacent shots and receivers. If these effects have not been corrected before 
noise removal then the noise removal process must presen/e them, one method for this is 
surface-cmsistent f-x prediction (Wang. 1996). 

Another applicatk)n of noise removal is common^ffeet or common-angle stacks for 
amplitude versus offset (AVO) or amplitude versus angle (AVA) analysis. Such stacks are used 
for the automatic computation of panimeters for interpretatton. These stacks require a low level 
of noise so the computed parameters are as accurate as possible. In 2-D acquisition, AVO/AVA 
stadcs form a 3-0 volume in whtoh the two spatial dimensions are CMP and either offset or 
angle. In the offset/angle dinoension there may be only one or two dozen traces, such that much 
of the data is on or near a spatial boundary. Oisadvantageously even the better noise removal 
methods, such as f-xy prediction, do not perf<mn well near spatial boundaries, r^uWng in 
distortion of the s^nal, and possible distortion of the computed parameters - creating the need 
for a ndse removal method ttiat performs well at spatial boundaries. 

Some conventional noise removal methods sudi as Karhunen-lx)eve. time-domain 
eigenimage filtering, Radon transform, and f-K and f-kk transforms, only work welt on plains- 
type data received from geologkial formattons in whteh nrast of the reflectors are flat 
Oisadvantageously m more structured data received from geotogtcal formations in which 
reflectors are strongly dipping or stoping, these methods become computationally expensh^e, 
difficult to use, or less effective. Accordingly, there is a need for a method of noise removal that 
perfonns well on ait types of geok)gy. 

Many methods for removing noise from seismic data are implemented using matrix 
operattons. For example, matrix compression is one process for determining an approximation 
to a given matrix, which approxinr)atk)n consumes less space than the original matrix. However, 
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matrix rank reduction is to determine the nearest (with respect to a particular matrix nonrn) 
Tsrk-K matrix to a given matrix. A nr^atrix norm measures the size of a matrix, for example, the 
"Frobenius norm' is the square root of the sum of the square of the matrix elements, whereas 
the ''LI matrix norm* is the sum of the absolute values of the elements of the given matrix. 

Prior art reviewed includes US 5,379,268 to Hutson, who teadies compression of the 
subject matrix (see daim 1 (b)) as the core of his improvement Compression and rank reduction 
have some similarities, but are not Menticdl, in that rank reduction can be used as a step within 
matrix compression, but rank reduction can also be performed without any resulting matrix 
compression. Hutson in 268 teaches compression - by active decomposition (expanding the 
original data matrix), then activeiy and selectively zeroing (whereas nrKXiifying to a value other 
than zero does not "compress") singular values in one of the resufb'ng matrices, thereafter 
recomposing those matrices into a single matrix that approximates the original - after which 288 
proceeds to teach (see daim 1(c)(1)) further processing of ttie signals approximated by the 
oompressed matrix. However, 268 is vague even about the kind of compression and processing 
perfomoed. For example modifying without zeroing cannot compress while zeroing inappropriate 
selections can actually be counterproductive by eliminating s^nal rather than noise. 

Oisadvantageously, existing noise removal algorithms do not handle erratic noise well 
For example, both SVO and Lanczos methods of matrix rank reduction attempt to find a rank-K 
matrix that is nearest to tr^ input matrix in a Frobenius-nomi sense. This is appropriate for 
removing random noise that has a Gaussian, or bell^haped, statistical distribution - which does 
not perform as well yAxen erratic non^Gaussian noise bursts are present. 

Processing seismic data is typicafly time consuming and expensive because it involves 
large quantities of data. Therefore, it is desirable to provide a solution to at least some of the 
above-described problems of the prior art redudng eitiier or both the quantity of data processed 
or the amount of processing required in relation to thai data. The prior art indudes some matrix 
rank reduction and frequency domain based metif^cds, but there is a need for a new combination 
of such metiiods that overcomes the above disadvantages, particulariy for ti\e purpose of noise 
suppresston in sections of seismic data. 
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SUMMARY OF THE INVENTION 

In accordance with one of its aspects the present invention comprises a method based 
on finding tfie rank K matrix nearest to a given matrix after which the smallest K value is used 
based on which the drfference plot show insignificant signs of signal. According to the method 
aspect of the present invention random noise is removed from a seismic data set by the 
following steps: the subject seismic data set is spatially divided into overlapping, redangular 
grids of traces; each said grid is processed by transforming it into the frequency domain using a 
Discrete Fourier Transfbmi (OFT); the grid is then separated into constant-frequency slices; all 
or a subset of said slices in matrix fomn are then individually rank-reduced after which the grid is 
refomied using the rank-reduced fomns; an inverse OFT is then executed on each (frequency 
ordered) trace of ttie reformed grid transfomning the traces back into the time domain but with 
noise suppressed and ready for use in vizualisation or interpretation; once all of the required 
grids are so processed, the entire seismic data set may be refomned using the refont\ed grids. 

in accordance with one of its aspects the present invention comprises a method of 
removing noise from a time*domdtn based CUBE of seismic data consisting of a plurality of 
Traces, the method comprising the steps: transform each said Trace into the frequency-domain, 
for the purpose of creating a frequency-domain based CUBE of seismic data, wherein the 
seismic data elements of said frequency-domain based CUBE are complex-valued; 
disassemble said frequency-domain based CUBE into a plurality of constant frequency slices, 
each of said constant frequency slices consisting of a plurality of seismic data elements; and for 
each constant frequenQr slice of said plurality of constant frequency slices: fomi one Matrix A m x 
n from each said constant frequency slice using said plurality of seismic data elements as the 
complex-valued elements of said Matrix A m x n ; rank-reduce Matrix A m x n to create a rank- 
reduced Matrix B m X n that is representative of Matrix A „ « «; substitute Matrix B ^ x n in place of 
Matrix A m I n . for the purpose of forming a proxy slice that is representative of sakJ constant 
frequency slice; assemble a proxy frequency-domain based cube using said proxy slice, for the 
purpose of accessing each proxy trace in a plurality of frequency ordered proxy traces that are 
representative of said plurality of Traces; and inverse transform into the time-don^in each proxy 
trace of said plurality of frequency ordered proxy traces, for the purpose of forming a noise- 
suppressed time*domain based proxy cube representative of said time-domain based CUBE of 
seismic data. 
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Processing seismic data via the method of the present invention results in a grid that is 
representative of an original grid, but advantageously has a better signal to noise ratio. This 
advantage results in part from having fewer non-trivial data elements in each of a set of 
representative grids, since me surviving data elements represent (he bulk of the composite 
signal related to genuine reflectors - whereas the trivial elements replace a large portion of the 
composite signal related to random noise. Unlike the prior art. there is no compression of the 
elements of the representative matrices. The method asped of the present invention improves 
the interpretability of processed seismic data and worics at almost any stage of signal 
processing. For example, since it preserves surface-consistent effects, this method may be 
applied before statics correction or deconvolutton. It may also be used to rennove ndse on 3-0 
volumes of stacked traces as well as comrrron-ofiset or common-angle stadcs. This method can 
be used to remove coherent noise from seismic traces and addresses non-unifomn shooting 
patterns. For stacked 3-0 volumes this method can be executed fiaster than f-xy pred{ctk)n 
filtering and advantageously provides better results along the boundary of Vtre subject volume. 
For removing noise on common-offset or comnm-angle stacks* and for removing noise on pre- 
stack data, the present invention is superior to standard f-xy prediction filtering since due to a 
synergy resulting from the manner in which the grid is formed, the method of the present 
invention is faster and can preserve surface-consistent ef^cHs hereby allowing it to be applied 
before statk^s con^on or deconvolution. For example, the ntethod of ttie present inventk>n 
uses an m-by-n 2-0 grid of seismic traces, the spatial locations of which need not be truly 
rectangular (e.g. a parallelc^nam). Si'milariy, the distances between grid lines in eittier the row or 
column directions need not be evenly spaced. And, the steps of the method of the present 
inventton need not be performed on an entire data set since an appropriate subset of ail 
available frequencies may be used to target reflectors located at specific co-ordinate pairs. The 
amount of noise removed using the method of the present invention can be increased by 
increasing the grid dimensions m and n. or by decreasing the rank K. A 2-D grid of seismic 
traces suitable for processing via the method of ttie present invention may. for example but not 
in limitation, originate from: 

a rectangle of traces extracted from a stacked 3-D volume. The trace grid being 

comprised of inline CMP bins in the row direction, and crosstine CMP bins in the column 

direction; 

traces from an unstatiked 2-D line. The grid is composed of comnion source 
traces in the row direction, and common receiver traces in the column direction; 
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traces from an unstacked 3-D volume, where the traces are taken from a single 
shot line and receiver line. The trace gnd being comprised of common source traces in 
the row direction, and common receiver traces in the column direc^n; or 

traces from common-offset or common-angle stacks for a sequence of CMPs. 
The trace grid being comprised of common-offeet or -angle traces in the row direction, 
and CMP traces in the column direction. 

In order to overconte the disadvantages of the prior art in one of its aspects the present 
invention comprises a novel method based on matrix rank reduction for removing noise from 
seismic data sets, whidi method is frequency domain based and less time consuming and 
expensive to execute by redudng the amount of data processing required. According to one 
embodiment of the method of the invention it is not necessary to fuJly decompose the subject 
matrix or to "zero ouT elements because partial decomposition results in the desired matrix 
elements being extracted early in the decomposition process, permitting decomposition to be 
terminated before completion. 

The accompanying drawings, which are incorporated in and constitute a part of this 
specification, illustrate preferred embodiments of the method, system, and apparatus according 
to the invention and, together with the description, serve to explain tiie principles of tiie 
invention. 

BRIEF DESCRIPTION OF THE DRAWINGS 

The present invention, in order to be easily understood and practised, is set out in the 
following non-limiting examples shown in the accompanying drawings, in which: 

FIG 1. is Prior Art and an illustration of the typical system used in the dcquisition of a 
single seismic shot 

FIG 2. is Prior Art and an illustration of a typical 2-0 stacked seismic section 

FIG 3. is Prior Art and an illustration of a typical 3-0 stacked section of seismic data the 
X and Y planar surface of which includes a grid of co-ordinate pairs the data associated 
with which grid over time fonrts a CUBE of seismic data 

FIG 4. is a summary fk>w chart of an embodiment of the method of the present invention 

FIG 5. is a more detailed flow chart illustrating the removal of noise from data associated 
with a grid of traces 
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F(G 10. i$ a more detailed flow chart illustrating one embodiment of the method of the 
present invention for the removal of noi$e fix>m a CUBE of seismic data 

FIG 6. is a flow chart itlustrating one way to rediice the rank of a matrix 

FIG 7. is a flow chart relating to the selection of rank K 

FIG 8. is a surface stacking diagram describing the acquisition of a 2-D acquisition 
seismic line 

FIG 9. illustrates the posittons of shots and receivers in a 3-0 acQuisitron seismic an^y, 
and the seiectton of a particular shot and receiver line for individual noise removal 

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT 

Reference is to be had to Figures 4 - 10 in which klenticai reference numbers Mentrfy 
Similar components* 

According to one embodiment of the method of the present invention Figure 4 outlines a 
method for the removal of random noise from a section of seismic data 41 0. The sectk>n may be 
in 2-D or 3-D. stacked or not stacked. Although there are no known theorettcal limits to the size 
of ttie grid that wouM benefit from sudi noise suppression, experience suggests that efficient 
results are achieved when grid sizes of between 20 x 20 and 30 x 30 traces are used. The 
reason for the difference betvveen theory and practice relates to the liade-off between noise 
suppresston and ^'mixing'* or "geotogy smearing**. For example, using any industry standard 
spadng and a section based on a grid of 1000 x 1000 traces, it would be benefidal to spatially 
divkle the larger data set into overiapping rectangular grids of 25 x 25 traoes. which grids 
overiap on all edges suc^ that the data set of the 1000^ section could be covered by 
approxin)dtety 2500 smaller grids each suitable for use in dusters of better focused subsets of 
grids targeting suspected reflectors of interest. While running the larger 1000^ data set as a 
single grid is possible and the results would indude less random noise, experience indicates 
that such results would be degraded by smearing or mixing the subsurface geology 
unnecessarily because reflectors are typically small enough that refieded energy from them is 
not relevant to the grid that may be created from an entire data set. Consequently, spatially 
dividing a data set into smaller grids covering a target of interest represents a balance between 
noise suppression and smearing. As the confidence level assodated with and interest in a 
particular target increases, according to an alternate embodiment of the method of the present 
invention, several smaller overiapping grids may be used to blanket the relevant taiget 
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coordinates resulting in multiple layers of overlap surrounding a central location of interest 
Experience shows that several smal)« dosety-assodated grids provide good conrelation and are 
well suited for noise suppression by matrix rank reduction. Therefore, as appropriate, section 
410 is spatially divided 420 into overlapping rectangular grids of traces. If a grid is missing 
traces (usually because it is located near the spatial boundary of said section), one inserts 
artificial traces having values of zero to complete it At step 430 noise is removed from each of 
the grids into which section 410 has been divided, which removal may be done in parallel or 
sequentially. At step 440 section 410 is **reformed" using the noise suppressed grids, which 
results in a section that is representative of the original. At step 450 the noise suppressed 
section may be used to interpret subsurface geology in place of section 410. 

Referring to Rgure 5 there is illustrated one embodiment of the sub-steps of step 430 
necessary to process each CUBE of data associated with the grids into which section 410 has 
been divided at step 420. Step 510 is the selection of a particular grid for noise suppression 
processing. At step 520 the transform of each trace In said grid ie taken in order to create a 
CUBE of seismic data in the frequency domain, which CUBE is then disassembled into constant 
frequency slices at step 530. As shown in steps 540 end 550 the resulting slices may be 
processed in parallel. Noise redu^on occurs when the data elements of each slice are placed 
into indhMual matrices that are subsequently rank-reduced. At step 560 a proxy cube 
n^presentative of the CUBE of the original grid is formed by using the rank-reduced matrices in 
place of the matrices formed at step 540. At step 570 the inverse transfonn is taken of ea^ 
trace in said pro)^ cube to return the ou4)ut to the time domain and create a noise suppressed 
proxy cube or "grid of traces* that may be used to interpret subsurface geology in place of the 
particular grid selected at step 510. 

Referring to Figure 10 there is Illustrated one embodiment of the method of the present 
invenfon according to which each selected "grid of Traces" (or CUBE, each for example having 
625 Traces) is input at step 51 1 for processing by first transforming at step 521 each Trace in 
the selected CUBE into the frequency domain using a Discrete Fourier Transfonn (OFT), 
thereby creating a record • across a finite range of frequencies - of the composite (>.a reflectron 
^ noise) signal sampled at each co-ordinate pair of the subject grid, in tiiis example ttie CUBE 
of 625 time traces becomes a transformed CUBE of 625 "traces** expressed in temns of 
frequen(y. such that the elements of the transfomied data set are complex-valued. An alternate 
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way of desoribing the shift between the time and frequency domains is the ''re-ordering'' of the 
grid ou^ut from time trace ordered to frequency ordered. 



At step 531 the frequency transfomned CUBE i$ sfroed rnto a plurality Qn this example 
"h**) of constant frequency slices. A slice is then taken out of the transforn^ CUBE 
representing the composite signal captured at each co-ordinate pair, but at a single frequency 
rather than time. In our example, each such 25 x 25 co-ordinate pair grid slice will include 625 
complex-valued data elements that are susceptible to processing via matrix math. For each 
slice(h) there is a oomplex-vatued Matrix A m x ^ of the same dimensions, such that at step 541 
we fomi one Matrix A(h) using the complex-valued data elements from tiie oonresponding 
sfice(h). As shown at step 570 the steps 541, 545, 551, 555, and 561 are repeated for eadi 
stice(h) corresponding to a selected value of h that may represent full or partial decomposition 
according to how much noise is to be removed from the grid of Traces or CUBE of interest 



According to a preferred embodiment of the method of the present invention, at step 545 
the resulting Matrix A(h) is decomposed by applying, for example. Singular Value 
Decomposition fSVD**) to expand it into left singular vectors, singular values, and right singular 
vectors - resulting in 3 separate matrices, one of which is en ordered diagonal matrix 
Conveniently, ttiese data elements result in a Matrix A, where: 

A mxn ^ U £ V^* where £Js an ordered diagonal matrix and U and V" are both unitary, 
since it is a property of SVO that U and are both unitary and that the elen^ents of I 
will be organized wim the largest at an. aaa. asa. . . . decreasing to 
It is a property of the seismic record creation process that true reflected signal tends to 
be recorded or present primarily in the higher vati^ elements (am a22» 333) of L 
whereas random noise signal tends to be dispersed along the entire diagonal of Z 

According to an ^emate embodiment of the method of ttie present invention step 545 
may be skipped and processing may proceed direcHy to rank-reduction of Matrix A(h} m x a. by 
alternate means for the purpose of aeating a rank-reduced complex-valued matrix of the same 
dimensions being Matrix B(h) that is representa^e of Matrix A(h) „ . 



According to a preferred embodiment of the method of the present invention, at step 551 
for each Matrix A(h) mxn. create a rank-reduced comptex-valued matrix of the same dimensions 
being Matrix B(h) of rank K. (istote: ea<^ Matrix B has fewer non-trivial elements than the Matrix 
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A that it represents) where K is less than the lesser of m or n. The selection of K is addressed 
t>elow. Matrix B(h) of rank K results when in the ordered diagonal matrix Z (resulting fronr) step 
620) ail but the top K elements along the diagonal are zeroed (as at step 630) and the matrix is 
refofmed (at step 640) as: B m * « * U JT V" . where T is the rank-reduced version of E having 
only the top K elements remaining non-trivial, it is an important distinction that the top K 
elements are not compressed or othenMse modified. Similariy, neither is the ordered d^nal 
matrix 27 compressed, but merely rank-reduced. It is to be understood that the purpose of this 
exercise is the suppression of noise in the output rather than the compresston of the data set of 
the CUBE to save storage space. 

According to an alternate embodiment it is contemplated that the top K elements may be 
weighted or otherwise adjusted or processed to advantage, p^tculariy where the related data 
elements an^ at the boundary of the subject CUBE from which they ere derived - or where there 
is an overlap of the coH^rdinate p^rs of such boundary with other CUBEs that relate to the 
target reflector(s) of interest. 

According to the present example, if the CUBE is sliced at intervals of 0.5 Hz between 0 
Hz and 250 Hz, processing wouM involve approximately 500 constant frequency slices each 
comprised of 625 complex-valued data elements. With 500 slices having been taken and 
processed. 500 rank reduced matrices B mxn - U ^ V" ^$uit and each may be used at step 
555 to fomi a frequency ordered proxy slice by substituting each Matrix B{h) for the 
corresponding Matrix A(h). 

At step 561 a frequency (domain based) ordered pro)^ cube Is formed using in this 
example each of the 500 rank reduced matrices B^Kn-Urv^dsa data set inserted at the 
frequency of the slice(h) that it represents. Once all 500 of the sltoes are so inserted the proxy 
cube will contain data lespecfing the composite signal recorded at each of the 625 co-ordinate 
pairs across a finite range of frequencies, however the ratio of true reflected to signal to noise 
signal will be increased. When h reaches a value that achteves the selected degree of 
decompositton (i.e. h s z) at step 570 the method proceeds to the final step 571 of the present 
embodiment and the frequency (domain based) ordered proxy cube is transfonned bade to the 
time domain. At step 571 an inverse DFT is perfbrmed on each of the 625 frequency ordered 
traces to return the data set to the time domain for visualization and other purposes. As each of 
these invent transforms are xecuted a proxy time trace representative of the corresponding 
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Original time trace associated with the subject co-ordinate pair is formed. Once all 625 proxy 
time traces are formed the proxy cut>e is available for correlating with overlapping and nearby 
proxy cubes that may be compared, summed, or othenvise processed for use interpreting the 
subsurface reflectors respecting which the original CUBEs were selected. 

Onoe all CUBEs are so noise-suppressed, one may fomi a representation of the entire 
original section by using all the proxy cubes. Where the co-ordinate pairs of the traces in neait>y 
grids are the same, such overlap zones may be addressed by summing grid traces at the same 
(co-ordinate pair) position after scaling them with "Vvelghts" so that the sum of the weights at the 
overiap position is one. A person of sktil in the art of seismic processing would know to select 
such weights to taper smoothly near the boundary of each rectangular grid. 

According to the method of the present invention the selectton of K for the purposes of 
rank reduction may be more refined or less discriminating (i*e. fine, moderate, coarse). A 
somewhat arbitrary choice of K 2 may for example to a moderate degree of noise 
suppression with a moderate degree of undesired toss of the signal ftiom genuine reflectors. 
Choosing K = 1 means that all but the an element of the ordered diagonal matrix £ wouM be 
zeroed leading to a harsh degree of noise suppresston, but almost certainly a significant and 
undesired loss of s^nai from genuine reflectors. By contrast, choosing K = 3 tends to leave 
genuine signal in tact while achieving less noise suppresston advantage. 

According to an attemate embodiment of the method of the present invention non- 
integer values of K may be applied to fine tune the signal to noise ratio of the result For 
example, a K value of 27 may be implemented by zeroing out ail but the ttiree largest singular 
values, and multiptying the third largest singular value by 0.7 before the matrix is recomposed. 
In this circumstance K no longer represents rank, but rather a degree of noise removal that is 
intemiediate that of rank 2 and rank 3. 

Since noise and signal are not universally distributed across all CUBEs it is 
advantageous to select K on a grid by grid basis penmitting a nr>ore refined processing with the 
goal of an optimal balance between noise suppression and signal dlstortkm. 

Figure 7 illustrates one way to select rank K to remove noise without distorting the 
reflected signal of a seismic data set For example for each of K values 1 through 5 (defined at 
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step 720) at step 740 calculate the difference between Matrix A (I.e. the input data) and the 
Matrix 6 of rank IC At step 750 plot the difference (to quickly visualize how much of the signal 
has t>een removed) between Matrix A and Matrix B for eadh value of K (i.e. 1 * 5) and at step 
760 select the value of K (.e.g. 3) the plot for which shows insignificant indicatkms of signal (i.e. 
looks random with IttUe coherence). A person of skill in the art of seismic processing will readily 
recognize viAm a difference plot contains too much signal 

Referring to Figure 6 there is illustrated one means for executing step 545 to decompose 
each Matrix A(h). A person of skill in mathematk:s would understand that SVD is only one 
to arrive at the rank reduced Matrix B m x « - U r , but that the dimenstons of the matrix are 
not to be changed by the mathematical process selected and the top K elen^nts maintain their ^ 
relationships. The inventor has proven that eigenimages may be used and summed more 
quickly, but tend to be less accurate than the results from using SVD for rank reducing Matrix 6 
mxn = U 27 V*^- For example, although the SVD means of fully decomposing a matrix wortcs well 
for rank redudion. reasonable approximattons to full decomposition can be computed using 
Lanczos bi-diagonatizatfon that can require as HtUe as one-tenth the computational time of the 
SVD method yet me quality of Lanczos results is similar to the SVD results. Advantageously, 
when removing noise from large data sets, the method of the present invention can be executed 
much (aster than the closest known competing method, being f-xy prediction filtering. 

The method of the present invention wortcs equally well on both flat and structured data 
because this meUiod does nothing to a noiseless seismic grid containing no more than K dips, 
which (unlike eigenimage flRering in the time domain) is because the method of the present 
invention operates in the frequency domain of each trace. 

According to a preferred embodiment of the method of the present inventbn not all of 
the frequency slices need to be rank reduced. Typically seismic traces are sampled in time at a 
rate such that the signal frequencies are a fraction of the Nyquist frequency. For example, it is 
common for seismic data to have significant signal only between frequencies 10 and 80 Hz (the 
appropriate signal band is well known to persons of skill in the art of seismic processing), yet the 
Nyquist frequency is often 125 or 250 Hz - consequently only matrices 540 based on frequency 
slk»s between 10 and 80 Hz need to be rank reduced. The remainder can be left unchanged or 
zeroed each resulting in a considerable savings in oomputatioa 
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Advantageously, according to an alternate embodiment, the method of the present 
invention can handle erratic noise quite well by identifying the rank K matrices that are near the 
subject input matrix in an Ll-nomi sense using a robust SVD (Hawkins. Uu, and Young,200l). 

As illustrated in the "surface stacking diagram" of Figure 8, the method of the present 
inventton can be applied to unstacked 2-D seismic data sets when unstacked traces 810 are 
law out on a two^imenstonal grid on which the trace shot (ordered by increasing 
receiver/station position) fbmis one axis 820 and the trace receiver fomis the other axis 830, 
such that the data has the appearance of a stacked 3-D section pemiitting noise removal to be 
perfonred as set out above. Referring to Figure 9 tt is illustrated that the method of the present 
inventton can be applied to an unstacked 3-0 seismic data set In a typical 3-0 acquisition, 
shots 910 are positioned spatially along a multitude of "shot lines", and receivers 920 are 
positioned spatially akmg a multitude of "receiver lines". To petfom noise removal according to 
an alternate application <A the method of the present inventton. extract all traces having been 
acquired on a single shot line 930 and a receiver line 940. These traces are then laM out on a 
spatial grid where shots from the shot line form one axis and receivers from the receiver line 
fomi the other axis - giving the data the appearance of a stacked 3-0 sectton on wNch noise 
removal may be performed as set out above. The foregoing process is repeated for all 
remaining combinattons of shot lines and receiver lines. 

The method of ttie present inventton works well for 2-0 and 3-D unstacked data sets 
because: the method is independent of x- and y-oonsistent statics (i,e. the Statics Property); the 
method is exact for a noiseless seismic grid that has no more than K dips, and has then had x- 
and y-conslstent filters applied (i.e. the Filtering Property); and if the method is exact Ibr a 
seismte grid then the method is also exact the same seismto grid whfeh has had rows or 
columns of traces duplicated or removed (i.e. the Shooting Properly). Advanlageouisly. as a 
result of the Stattos and FBtering properties, and tile fact that the nriatrw rows and columns are 
selected to represent common shots and receivers, the random noise removal can presen/e 
surface-consistent effects, altowing the method to be applied at a very eariy stage of 
processing. To extract r»:tangular grids from unstacked 2-0 and 3-0 data sets for noise 
removal, the x axis represents shots and y represents receivers because then surface- 
consistent (that is. shot and receiver) effects are left undistorted by the method as a result of a 
synergy between the method's ability to absort). or leave undistorted. x- and y-consistent 
effects, and the manner of extracting rectangular grids of traces from pre-stack data sets. 
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Advantageously, the method of the present invention works well along a straight spatial 
boundary, since from ttie method's point of view there is no boundary, which makes the method 
welkuited for removir^ ncise from common-offset or common*angle stacks, in whidi many of 
the traces are at or near a boundary. For comnrK^n-offset or common-angle stacks from a 2-D 
acquisition, the traces are naturally laid out in a 2-D spatial grid, making it possible to perform 
noise removal as if it were a stacked ^0 sectk)n. Advantageously, the n^thod is independent 
of the row and column ordering as a result of the Ordering Proper^. 

According to an alternate embodiment of the method of the present inventfon noise 
reduction can be designed on one set of data, but applied on another The design data can be 
taken from different time windows of the same traces as the applrcation data, or from a different 
set of traces. This is made possible where matrix A hokls the DFT values for a given frequency 
of the design data, and matrix C holds the DFT values for a given frequency of the application 
data * it is possible to calculate matrix B by projecting matrix C onto the rank K subspace of 
matrix A corresponding (o its first K singular values. 

According to an alternate embodiment of the method of the present invention, in the 
frequency x-y plane for a given frequency a rank K matrbc may be produced using eigen- 
analysis wherein K is the number of plane waves, which fad allows tirte separation of plane 
waves by eigen-image decomposition. A single frequency slice is rank-reduced by placing this 
2-D grid of complex DFT values into a complex-valued matrix of tiie same dimensions, finding 
the nearest rank-K matrix to this matrix, where K is some value greater than or equal to one. 
and replacing the constant-frequency slice values with the values from the rank-K matrix. 

According to an alternate embodiment of the mettK>d of the present invention, by 
applying different noise filtens to the design data, it is possible to remove coherent noise from 
seismic data, as well as random noise, whidi permits tattering the signal subspace to avoid, and 
thus renrK)ve, coherent energy. 

Although the disclosure describes and illustrates various embodiments of the invention, 
it is to be understood ttiat the invention is not limited to these particular embodiments. Many 
variations and modifications will now occur to those skilled in the art of processing seismic data. 
For full definition of the scope of the invention, reference is to be made to the appended daints. 



